Problem B: Short-term effect of air pollution on mortality

Complete pooling model

Motivated by the above relationships I start with the predictors \(SO2\) and \(TSP\):

\[ \text{avg_mort}_i \sim \mathcal{N}(\mu_i,\sigma^2)\\ \mu_i = \alpha + \beta\; \text{SO2}_i + \gamma \; TSP_i \] which has \(\mathcal{N}(0,30)\) prior distributions on \(\alpha,\beta,\gamma\) and \(Cauchy(0,30)\) on \(\sigma\).

## arrange data into a list
data <- list( 
  N = nrow(milan.mort), 
  avg_mort = milan.mort$avg.mort,
  SO2 = milan.mort$SO2+50,
  TSP = milan.mort$TSP+30,
  mean_temp = milan.mort$mean.temp
)

## compile the model
model1 <- stan_model("model1.stan", auto_write = TRUE)

## fit the model
model1_fit <- sampling(model1, data = data, cores=4)
saveRDS(model1_fit, "model1_fit.rds")

# Load fitted model
model1_fit <- readRDS("model1_fit.rds")

print(model1_fit, pars = c('alpha','beta','gamma','sigma'))

## Inference for Stan model: model1.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## alpha 28.17    0.01 0.19 27.78 28.04 28.17 28.29 28.53  1157    1
## beta   0.03    0.00 0.00  0.03  0.03  0.03  0.03  0.03  4235    1
## gamma -0.01    0.00 0.00 -0.01 -0.01 -0.01 -0.01  0.00  1569    1
## sigma  4.72    0.00 0.06  4.61  4.68  4.72  4.76  4.83  1564    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun 18 09:29:31 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

## traceplots
plot(model1_fit, plotfun = "trace", pars = c('alpha','beta','gamma','sigma')) + ggtitle("traceplots") 

## acf plots
mcmc_acf(as.matrix(model1_fit, pars=c("alpha","beta","gamma","sigma")))

# for model comparison
log_lik_1 <- extract_log_lik(model1_fit)
loo_1 <- loo(log_lik_1)

## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and 
## MCSE estimates will be over-optimistic.

Rhat value equal to one indicates the convergence of the chains, which is confirmed by the traceplots mixing well. Also, the autocorrelation decreases for all parameters.

# relationship bw estimated parameters
mcmc_scatter(as.matrix(model1_fit, pars = c('alpha','beta')), alpha = 0.2)

mcmc_scatter(as.matrix(model1_fit, pars = c('beta','gamma')), alpha = 0.2)

mcmc_scatter(as.matrix(model1_fit, pars = c('alpha','sigma')), alpha = 0.2)

The scatterplot shows a relationship between variables beta and gamma.

# real vs estimated checks
y_rep <- as.matrix(model1_fit, pars = "y_rep")
y <- milan.mort$avg.mort

pp_check_density(y_rep, y=y)

ppc_stat(y = y, yrep = y_rep, stat = "mean")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ppc_stat(y = y, yrep = y_rep, stat = "sd")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

In the posterior predictive check, the simulated distributions resemble the original one.

The residuals look mostly positive, this means that the model understimates the real outcomes.

Including temperature in the model

ggplot(milan.mort, aes(x = mean.temp, y = milan.mort$avg.mort)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)

The plot shows an evident relationship between the average number of deaths and the mean daily temperature, so I include the information about temperature in the model \[ \text{avg_mort}_i \sim \mathcal{N}(\mu_i,\sigma^2)\\ \mu_i = \alpha + \beta\; \text{SO2}_i + \gamma \; TSP_i + \delta \; \text{mean_temp}_i \] I’m using the same priors from the first model, and a \(\mathcal{N}(0,30)\) on \(\delta\).

## compile the model
model2 <- stan_model("model2.stan", auto_write = TRUE)

## fit the model
model2_fit <- sampling(model2, data = data, cores = 4)
saveRDS(model2_fit, "model2_fit.rds")

# Load fitted model
model2_fit <- readRDS("model2_fit.rds")

print(model2_fit, pars = c('alpha','beta','gamma','delta','sigma'))

## Inference for Stan model: model2.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## alpha 34.06    0.01 0.32 33.43 33.85 34.07 34.27 34.67  1103    1
## beta   0.02    0.00 0.00  0.02  0.02  0.02  0.02  0.02  2053    1
## gamma -0.01    0.00 0.00 -0.01 -0.01 -0.01 -0.01 -0.01  3155    1
## delta -0.27    0.00 0.01 -0.29 -0.28 -0.27 -0.26 -0.25   946    1
## sigma  4.43    0.00 0.05  4.33  4.39  4.43  4.46  4.54  1629    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun 18 09:32:32 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

# for model comparison
log_lik_2 <- extract_log_lik(model2_fit)
loo_2 <- loo(log_lik_2)

## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and 
## MCSE estimates will be over-optimistic.

## pp check
y_rep <- as.matrix(model2_fit, pars = "y_rep")
y <- milan.mort$avg.mort

pp_check_density(y_rep, y)

std_residuals(y_rep, y)

out_perc(y_rep, y)

## [1] 0.01398026

ppc_stat(y = y, yrep = y_rep, stat = "mean")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ppc_stat(y = y, yrep = y_rep, stat = "sd")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# grouped data
ppc_stat_grouped(
  y = y, 
  yrep = y_rep, 
  group = milan.mort$day.of.week, 
  stat = 'mean',
  binwidth = 0.2
)

ppc_stat_grouped(
  y = y, 
  yrep = y_rep, 
  group = milan.mort$day.of.week, 
  stat = 'sd',
  binwidth = 0.2
)

Even if the chains are converging correctly and the autocorrelation is decreasing, the predictions on the mean value are not always plausible over the different days of the week.

Hierarchical model

We could take into account the hierarchical structure of the data in week days and try modelling the variation across them. In this case the response variable is tot.mort, because it would make no sense to model daily variations on averages calculated along the last three days. Both the intercept and the standard deviation are estimated for each different day of the week, moreover this time the priors are \(\mathcal{N}(0,10)\).

\[ \text{avg_mort}_{d} \sim \mathcal{N}(\mu_{d},\sigma_d)\\ \mu_{d} = \alpha_d + \beta\; \text{SO2} + \gamma \; TSP + \delta \; \text{mean_temp} \]

# reload data
data(milan.mort)

# hierarchical data
hier_data <- list(
  N = nrow(milan.mort), 
  D = length(unique(milan.mort$day.of.week)),
  tot_mort = milan.mort$tot.mort,
  SO2 = milan.mort$SO2+50,
  TSP = milan.mort$TSP+30,
  mean_temp = milan.mort$mean.temp,
  day_of_week = milan.mort$day.of.week,
  holiday = milan.mort$holiday
)
str(hier_data)

## List of 8
##  $ N          : int 3652
##  $ D          : int 7
##  $ tot_mort   : num [1:3652] 45 32 37 33 36 45 46 38 29 39 ...
##  $ SO2        : num [1:3652] 317 425 326 490 404 ...
##  $ TSP        : num [1:3652] 140 183 192 228 265 ...
##  $ mean_temp  : num [1:3652] 5.6 4.1 4.6 2.9 2.2 0.7 -0.6 -0.5 0.2 1.7 ...
##  $ day_of_week: int [1:3652] 2 3 4 5 6 7 1 2 3 4 ...
##  $ holiday    : int [1:3652] 1 0 0 0 0 0 0 0 0 0 ...

## compile the model
model3 <- stan_model("model3.stan", "auto_write" = TRUE)

## fit the model
model3_fit <- sampling(model3, data = hier_data, iter = 1000, chains = 4, cores=2)
saveRDS(model3_fit, "model3_fit.rds")

# Load fitted model
model3_fit <- readRDS("model3_fit.rds")

print(model3_fit, pars = c('alpha','beta','gamma','delta','sigma'))

## Inference for Stan model: model3.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##           mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## alpha[1] 32.83    0.02 0.57 31.72 32.44 32.84 33.22 33.90   875    1
## alpha[2] 32.99    0.02 0.57 31.86 32.61 32.99 33.37 34.10   886    1
## alpha[3] 32.41    0.02 0.56 31.34 32.01 32.42 32.80 33.44  1082    1
## alpha[4] 32.69    0.02 0.58 31.56 32.31 32.69 33.08 33.80   777    1
## alpha[5] 32.39    0.02 0.58 31.23 32.00 32.39 32.79 33.53   916    1
## alpha[6] 32.16    0.02 0.55 31.08 31.78 32.17 32.55 33.22   787    1
## alpha[7] 32.17    0.02 0.52 31.15 31.82 32.17 32.53 33.19   868    1
## beta      0.02    0.00 0.00  0.02  0.02  0.02  0.02  0.02  1222    1
## gamma    -0.01    0.00 0.00 -0.01 -0.01 -0.01  0.00  0.00  2038    1
## delta    -0.23    0.00 0.02 -0.26 -0.24 -0.23 -0.21 -0.19   772    1
## sigma[1]  7.14    0.00 0.22  6.71  6.98  7.14  7.29  7.57  2112    1
## sigma[2]  6.72    0.00 0.21  6.33  6.57  6.71  6.86  7.17  1855    1
## sigma[3]  6.52    0.00 0.20  6.15  6.39  6.52  6.66  6.92  2179    1
## sigma[4]  6.85    0.00 0.21  6.47  6.71  6.84  6.99  7.26  2152    1
## sigma[5]  6.99    0.01 0.22  6.58  6.83  6.98  7.14  7.43  1921    1
## sigma[6]  6.59    0.00 0.20  6.21  6.45  6.59  6.73  7.00  2143    1
## sigma[7]  6.68    0.01 0.21  6.30  6.53  6.67  6.82  7.13  1389    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun 18 09:37:08 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

#computing psis-looic
log_lik_3 <- extract_log_lik(model3_fit)
loo_3 <- loo(log_lik_3)

## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and 
## MCSE estimates will be over-optimistic.

# areas plots
mcmc_areas(as.matrix(model3_fit, pars = c('alpha')))

mcmc_areas(as.matrix(model3_fit, pars = c('sigma')))

## traceplots
plot(model3_fit, plotfun = "trace", pars = c('alpha','beta','gamma','sigma')) + ggtitle("traceplots") 

## acf plots
mcmc_acf(as.matrix(model3_fit, pars=c("alpha","beta","gamma","sigma")))

## parameters distributions
mcmc_hist(as.matrix(model3_fit, pars = c('alpha','sigma')))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## pp checks
y_rep <- as.matrix(model3_fit, pars = "y_rep")
y <- milan.mort$tot.mort

pp_check_density(y_rep, y)

# real mean vs estimated mean
ppc_stat_grouped(
  y = y, 
  yrep = y_rep, 
  group = milan.mort$day.of.week, 
  stat = 'mean',
  binwidth = 0.2
)

ppc_stat_grouped(
  y = y, 
  yrep = y_rep, 
  group = milan.mort$day.of.week, 
  stat = 'sd',
  binwidth = 0.2
)

Now the standardized residuals look and the model is able to correctly describe the average number of deaths among different days of the week.

# add month information
milan.mort <- milan.mort %>%
  mutate(date = as.Date("1980-01-01") + 0:(length(milan.mort$day.num)-1)) %>% 
  mutate(month = month(date))

# hierarchical data
hier_data <- list(
  N = nrow(milan.mort), 
  D = length(unique(milan.mort$month)),
  tot_mort = milan.mort$tot.mort,
  SO2 = milan.mort$SO2+50,
  TSP = milan.mort$TSP+30,
  mean_temp = milan.mort$mean.temp,
  day_of_week = milan.mort$month,
  holiday = milan.mort$holiday
)
str(hier_data)

## List of 8
##  $ N          : int 3652
##  $ D          : int 12
##  $ tot_mort   : num [1:3652] 45 32 37 33 36 45 46 38 29 39 ...
##  $ SO2        : num [1:3652] 317 425 326 490 404 ...
##  $ TSP        : num [1:3652] 140 183 192 228 265 ...
##  $ mean_temp  : num [1:3652] 5.6 4.1 4.6 2.9 2.2 0.7 -0.6 -0.5 0.2 1.7 ...
##  $ day_of_week: num [1:3652] 1 1 1 1 1 1 1 1 1 1 ...
##  $ holiday    : int [1:3652] 1 0 0 0 0 0 0 0 0 0 ...

## fit the model
model3_fit2 <- sampling(model3, data = hier_data, iter = 1000, chains = 4, cores=4)
saveRDS(model3_fit2, "model3_fit2.rds")

# Load fitted model
model3_fit2 <- readRDS("model3_fit2.rds")

print(model3_fit2, pars = c('alpha','beta','gamma','delta','sigma'))

## Inference for Stan model: model3.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##            mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## alpha[1]  31.80    0.02 0.65 30.53 31.35 31.80 32.24 33.02  1764    1
## alpha[2]  31.45    0.02 0.67 30.13 31.00 31.45 31.92 32.76  1784    1
## alpha[3]  31.40    0.02 0.61 30.24 30.99 31.41 31.82 32.55  1155    1
## alpha[4]  28.46    0.02 0.62 27.28 28.06 28.46 28.87 29.68  1290    1
## alpha[5]  27.45    0.02 0.73 26.07 26.95 27.47 27.92 28.92  1270    1
## alpha[6]  26.02    0.02 0.86 24.43 25.42 26.04 26.59 27.72  1193    1
## alpha[7]  23.27    0.03 1.01 21.33 22.56 23.26 23.95 25.29   991    1
## alpha[8]  19.54    0.03 0.94 17.76 18.87 19.55 20.19 21.40  1054    1
## alpha[9]  22.90    0.02 0.81 21.37 22.35 22.88 23.45 24.54  1136    1
## alpha[10] 27.82    0.02 0.68 26.45 27.33 27.83 28.28 29.17  1145    1
## alpha[11] 28.68    0.02 0.60 27.54 28.27 28.69 29.10 29.85  1521    1
## alpha[12] 30.01    0.01 0.59 28.85 29.62 30.01 30.42 31.12  1845    1
## beta       0.02    0.00 0.00  0.02  0.02  0.02  0.02  0.02  2511    1
## gamma     -0.01    0.00 0.00 -0.01 -0.01 -0.01 -0.01  0.00  2539    1
## delta      0.16    0.00 0.04  0.08  0.13  0.16  0.18  0.23  1158    1
## sigma[1]   6.88    0.01 0.27  6.34  6.69  6.87  7.06  7.41  2325    1
## sigma[2]   8.02    0.01 0.34  7.37  7.80  8.01  8.26  8.71  2869    1
## sigma[3]   6.90    0.01 0.28  6.37  6.71  6.89  7.09  7.46  2764    1
## sigma[4]   5.91    0.00 0.24  5.46  5.74  5.91  6.07  6.42  2528    1
## sigma[5]   5.77    0.00 0.24  5.33  5.60  5.76  5.93  6.25  2604    1
## sigma[6]   5.99    0.01 0.24  5.55  5.82  5.98  6.15  6.48  1330    1
## sigma[7]   7.96    0.01 0.32  7.37  7.74  7.95  8.18  8.61  1963    1
## sigma[8]   6.31    0.01 0.25  5.86  6.14  6.30  6.48  6.84  2039    1
## sigma[9]   5.57    0.01 0.23  5.13  5.41  5.56  5.72  6.06  1987    1
## sigma[10]  6.19    0.01 0.25  5.71  6.02  6.19  6.35  6.71  1936    1
## sigma[11]  6.22    0.00 0.25  5.74  6.04  6.20  6.38  6.76  2693    1
## sigma[12]  6.09    0.00 0.25  5.63  5.92  6.08  6.24  6.61  2693    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Jun 18 09:46:42 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

#computing psis-looic
log_lik_4 <- extract_log_lik(model3_fit2)
loo_4 <- loo(log_lik_4)

## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and 
## MCSE estimates will be over-optimistic.

# areas plots
mcmc_areas(as.matrix(model3_fit2, pars = c('alpha')))

mcmc_areas(as.matrix(model3_fit2, pars = c('sigma')))

## traceplots
plot(model3_fit2, plotfun = "trace", pars = c('alpha','beta','gamma','sigma')) + ggtitle("traceplots") 

## acf plots
mcmc_acf(as.matrix(model3_fit2, pars=c("alpha","beta","gamma","sigma")))

## parameters distributions
mcmc_hist(as.matrix(model3_fit2, pars = c('alpha','sigma')))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## pp checks
y_rep <- as.matrix(model3_fit2, pars = "y_rep")
y <- milan.mort$tot.mort

pp_check_density(y_rep, y)

std_residuals(y_rep, y)

out_perc(y_rep, y)

## [1] 0.07667032

ppc_stat_grouped(
  y = y, 
  yrep = y_rep, 
  group = milan.mort$month, 
  stat = 'mean',
  binwidth = 0.2
)

ppc_stat_grouped(
  y = y, 
  yrep = y_rep, 
  group = milan.mort$month, 
  stat = 'sd',
  binwidth = 0.2
)